Effect of charged impurity correlation on transport in monolayer and bilayer graphene 
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We study both monolayer and bilayer graphene transport properties taking into account the 
presence ol correlations in the spatial distribution of charged impurities. In particular we find 
that the experimentally observed sublinear scaling of the graphene conductivity can be naturally 
explained as arising from impurity correlation effects in the Coulomb disorder, with no need to 
assume the presence of short-range scattering centers in addition to charged impurities. We find 
that also in bilayer graphene correlations among impurities induce a crossover of the scaling of the 
conductivity at higher carrier densities. We show that in the presence of correlation among charged 
impurities the conductivity depends nonlinearly on the impurity density m and can even increase 
with m. 



I. INTRODUCTION 

The scaling of the conductivity a as a function of gate- 
voltage, proportional to the average carrier density n, is 
invaluable in characterizing the properties of graphcnoi. 
The functional dependence of <r(n) at low tempera- 
tures contains information 2 ^ about the nature of disor- 
der in the graphene environment (i.e., quenched charged 
impurity centers, lattice defects^, interface roughness^, 
ripples 6 -^, resonant scattering centers^—, etc.) giving 
rise to the dominant scattering mechanism. At finite 
temperatures electron-phonon scattering contributes to 
the resistivity^ 2 ^—. However, in graphene the electron- 
phonon scattering is very weak and it becomes impor- 
tant only at relatively high temperatures (> 400 A'), 
as evidence also from the fact that around room tem- 
perature the temperature dependence of a appears to 
be dominated by activation processe s 15 ' 16 . The quan- 
titative weakness of the electron-phonon interaction in 
graphene gives particular impetus to a thorough under- 
standing of the disorder mechanisms limiting graphene 
conductivity since this may enable substantial enhance- 
ment of room temperature graphenc-based device for 
technological applications. This is in sharp contrast 
to other high-mobility 2D systems such as GaAs-based 
devices whose room-temperature mobility could be or- 
ders of magnitude lower than the corresponding low- 
temperature disorder-limited mobility due to strong car- 
rier scattering by phononsii. Therefore, a complete un- 
derstanding of the disorder mechanisms controlling o~(n) 
in graphene at T = is of utmost importance both from 
a fundamental and a technological prospective. 

The experimental study of er(n) in gated graphene 
goes back to the original discovery of 2D grapheneji^ 
and is a true landmark in the physics of electronic ma- 
terials. Essentially, all experimental work on graphene 
begins with a characterization of <j(n) and the mobil- 
ity, p = a/(ne). A great deal is therefore knownii^r— 
about the experimental properties of cr(n) in graphene. 
The most important features of the experimentally ob- 
served cr(n)i2r— in monolayer graphene (MLG) are: (1) a 



non-universal sample-dependent minimum conductivity 
o~(n f=s 0) = <T m i n at the charge neutrality point (CNP) 
where the average carrier density vanishes; (2) a linearly 
increasing, o~(n) oc n , conductivity with increasing car- 
rier density on both sides of the CNP upto some sample 
dependent characteristic carrier density; (3) a sublinear 
o~(n) for high carrier density, making it appear that the 
very high density cr(n) may be saturating. 

To explain the above features of o~(n) a model has been 
proposed 2 ^ - — with two distinct scattering mechanisms: 
the long-range Coulomb disorder due to random back- 
ground charged impurities and static zero-range (often 
called "short-range") disorder. The net graphene con- 
ductivity with these two scattering sources is then given 
by a = p _1 = (p c + p s ) , where p c and p s are resistivi- 
ties arising respectively from charged impurity and short- 
range disorder. It has been shown that 2 -! 2 ^— p c ~ l/n 
and p s ~ constant in graphene, leading to o~(n) going as 



where the density independent constants A and C are 
known 2 - as functions of disorder parameters; A, arising 
from Coulomb disorder, depends on the impurity density 
(rii) (and also weakly on their locations in space) and the 
background dielectric constant (k) whereas the constant 
C, arising from the short-range disorder^!, depends on 
the strength of the white-noise disorder characterizing 
the zero-range scattering. Eq. ([T]) clearly manifests the 
observed o~(n) behavior of graphene for n ^ since 
a(n < A/C) ~ n, and a(n > A/C) ~ 1/C with a(ri) 
showing sublinear (C + A/n)^ 1 behavior for n ~ A/C. 

The above-discussed scenario for disorder-limited 
graphene conductivity, with both long-range and short- 
range disorder playing important qualitative roles at in- 
termediate (rii < n ^ A/C) and high (n > A/C) carrier 
densities respectively, has been experimentally verified by 
several groups ^ 22 ' 24 . There is, however, one serious is- 
sue with this reasonable scenario: although the physical 
mechanism underlying the long-range disorder scattering 
is experimentally established 2 - ' 19 ' 20 to be the presence of 
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unintentional charged impurity centers in the graphene 
environment, the physical origin of the short-range dis- 
order scattering is unclear and has so far eluded direct 
imaging experiments. As a matter of fact the experimen- 
tal evidence suggests that point defects (e.g. vacancies) 
are rare in graphene and should produce negligible short- 
range disorder. There have also been occasional puzzling 
conductivity measurements [e.g., Ref. [30ll3l| reported in 
the literature which do not appear to be explained by the 
standard model of independent dual scattering by long- 
and short-range disorder playing equivalent roles. 

Recently a novel theoretical model has been 
proposed^ 2 - that is able to semiquantitativcly explain all 
the major features of cr(n) observed experimentally as- 
suming only the presence of charged impurities. The key 
insight on which the model relies is the fact that in ex- 
periments, in which the samples are prepared at room 
temperature and are often also current annealed, it is 
very likely that spatial correlations arc present among 
the charged impurities. In particular this model is able 
to explain the linear (sublinear) scaling of <r(n) in MLG 
at low (high) n without assuming the presence of short- 
range scattering centers. 

In this work we first review the transport model pro- 
posed in Ref. [Hj], an d then extend it to the case of 
bilayer graphene (BLG). We find that, as in MLG, the 
presence of spatial-correlations among impurities is able 
to explain a crossover of the scaling of cr(n) from low 
n to high n in BLG, as observed in experiments, and 
that, because of the spatial correlations, a depends non- 
monotonically on the impurity density Hi. 

The remainder of this paper is structured as follows. 
In Section |TT] we present the model and the results for the 
structure factor S*(q) that characterizes the impurity cor- 
relations. With the structure factor calculated in Sec. |TT] 
we provide the transport theory in Section IIIII and Sec- 
tion IIVI In Section IIIII wc study the density-dependent 
conductivity u(n) of monolayer graphene in the presence 
of correlated charged impurities. We calculate a(n) at 
higher carrier density using the Boltzmann transport the- 
ory. We also evaluate cr(n) applying both Thomas-Fermi- 
Dirac theory^ and effective medium theory^ to charac- 
terize the strong carrier density inhomogeneities close to 
the charge neutrality point. In Section [IV! we apply the 
Boltzmann transport theory and the effective medium 
theory for correlated disorder to bilayer graphene and 
discuss the qualitative similarities and the quantitative 
differences between monolayer and bilayer graphene. We 
briefly review the experimental situation in Section [V] 
We then conclude in Section [VTI 



II. STRUCTURE FACTOR S(q) OF 
CORRELATED DISORDER 

In this section we describe the model used to calculate 
the structure factor S(q) for the charged impurities. We 
then present results for S*(q) obtained using this model 



via Monte Carlo simulations. The Monte Carlo results 
are then used to build a simple continuum approximation 
for S'(q), which captures all the features of S*(q) that are 
relevant for the calculation of a(n). 



A. Model for the structure factor S(q) 

To calculate S(q) we follow the procedure presented in 
Ref. HH, adapted to the case of a honeycomb structure. 
The approach was applied to study the effects of impu- 
rity scattering in GaAs hetero junctions and successfully 
explained the experimental observation of high- mobilities 
(e.g. greater than 10 7 cm 2 /(V-s)) in modulation-doped 
GaAs hcterostructures. The possible charged impurity 
positions on graphene form a triangular lattice speci- 
fied by ylm = + hAI. The vectors a = (l,0)ao 
and b = (^3/2,1/2)00 defined in the x-y plane, with 
ao = 4.92A, which is two times the graphene lattice con- 
stant since the most densely packed phase of impurity 
atoms (e.g. K as in Ref. |20| ) on graphene is likely to 
be an m x m phase with m = 2 for K— . The structure 
factor, including the Bragg scattering term, is given by 
the following equation: 

^^(E 6 " 1 '^) ( 2 ) 

where , Tj are the random positions on the lattice i"lm 
of the charged impurities and the angle brackets denote 
averages over disorder realizations. Introducing the frac- 
tional occupation / = Ni /N of the total number of avail- 
able lattice sites N by the number of charged impurities 
Ni, and the site occupation factor e^M equal to 1 if site 
r/ is occupied or zero if unoccupied, we can rewrite Eq. 
© as 

J LM 

in which the sum is now over all the available lattice 
sites (not only the ones occupied by the impurities). By 
letting Clm = {(-lai^o) I P we can rewrite Eq. (|3|) as: 

s(q) = /E c ^ e "™- ( 4 ) 

LM 

We then subtract the Bragg scattering term from this 
expression considering that it does not contribute to the 
resistivity obtaining 

5(q) = /^(C iM -l)e iq - r ^. (5) 

LM 

It is straightforward to see that for the totally random 
case, the structure factor is given by S(q) = 1 — / and 
n.i ~ 4.8/ x 10 14 cm -2 . For the correlated case we assume 
that two impurities cannot be closer than a given length 



FIG. 1: (a) Density plot of the structure factor S(q) obtained 
from Monte Carlo simulations for oo = 4.92 A and ro = 5ao. 
(a) m = 0.95 x 10 12 cm" 2 ; (b) m = 4.8 x 10 12 cm" 2 . 



fo < ft = (Tr^i)" 1 ^ 2 defined as the correlation length. 
This model is motivated by the fact that two charged im- 
purities cannot be arbitrarily close to each other because 
the Coulomb repulsion among the impurities during de- 
vice growth and there must be a minimum separation 
between them. 



B. Monte Carlo results for 5(q) 



FIG. 2: (a) and (b) show the calculated structure factor S(q) 
for two values of impurity density rn. (a) rn = 0.95 x 10 12 
cm" 2 ; (b) Hi = 4.8 x 10 12 cm" 2 . The solid lines show S(q) 
using Eq. (JHJ . Dot-dashed and dashed lines show the Monte 
Carlo results for two different directions of q from x-axis, 
9 — and 9 — 30°, respectively. 



introduce a simple continuum model being able to repro- 
duce for small q the structure factor obtained via Monte 
Carlo simulations. A reasonable continuum approxima- 
tion to the above discrete lattice model is given by the 
following pair distribution function g(r) (r is a 2D vector 
in the graphene plane), 



Using Monte Carlo simulations carried out on a 200 x 
200 triangular lattice with 10 6 averaging runs and peri- 
odic boundary conditions we have calculated the struc- 
ture factor given by Eq. (JSJ. In the Monte Carlo cal- 
culation a lattice site is chosen randomly and becomes 
occupied only if it is initially unoccupied and has no 
nearest neighbors within the correlation length ro- This 
process is repeated until the required fractional occupa- 
tion for a given impurity density is obtained. Once the 
configuration is generated, the Clm can be numerically 
determined after doing the ensemble average. In the nu- 
merical calculations, we use only statistically significant 
Clm, i-e., \ylm — r oo| < 3r , since Clm is essential unity 
for \y L m - rool > 3r . 

In Fig. [TJ we present a contour plot of the structure 
factor S(q) obtained from the Monte Carlo simulations 
for two different values of the impurity density. For ro 7^ 
the structure factor is suppressed at small momenta. 
Moreover the suppression of S'(q) at small momenta is 
more pronounced, for fixed ro, as rii is increased as it 
can be seen comparing the two panels of Fig. [1] The 
magnitude of 5(q) at small q mostly determines the d.c. 
conductivity and therefore, from the results of Fig. [TJ is 
evident that the presence of spatial correlations among 
the charged impurities will strongly affect the value of 
the conductivity. 



ff( r ) 



|r| < r 

1 |r| > r 



(6) 



for the impurity density distribution. In terms of the pair 
correlation function g(r) the structure factor is given by: 



S(q) = 1 



"[3(r)-l] 



(7) 



For uncorrelated random impurity scattering, as in the 
standard theory, g(r) = 1 always, and S(q) = 1. With 
Eqs. © and 0, we have 



S(q) 



1 



2nn i — J 1 (qro) 

q 



(8) 



where J\(x) is the Bessel function of the first kind. Fig. 
[5] shows 5(q) obtained both via Monte Carlo simula- 
tions and by using the simple continuum analytic model 
[Eq. ((SJ] for a few values of r and n.;. We can see that 
the continuum model reproduces extremely well the de- 
pendence of the structure factor on q for small momenta, 
i.e. the region in momentum space that is relevant for 
the calculation of a. 



III. MONOLAYER GRAPHENE 
CONDUCTIVITY 



C. Continuum model for S(q) 

Given that the value of the d.c. conductivity depends 
almost entirely on the value of ^(q) at small momenta, 
as discussed in Sections Mil and IIV1 it is convenient to 



In this section, we explore how the spatial correlations 
among charged impurities affect monolayer graphene 
transport properties. To minimize the parameters en- 
tering the model we assume the charged impurities to be 
in a 2D plane placed at an effective distance d from the 
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graphcne sheet (and parallel to it). 

We first study the density-dependent conductivity in 
monolayer graphene transport for large carrier densi- 
ties (n 3> ni) using the Boltzmann transport theory, 
where the density fluctuations of the system can be ig- 
nored. We then discuss cr(ro) close to the charge neutral- 
ity point, where the graphene landscape breaks up into 
puddle a 33 i 36 ~ — of electrons and holes due to the effect of 
the charged impurities using the effective medium theory 
developed in Ref . pt| . 



A. High density: Boltzmann transport theory 

Using the Boltzmann theory for the carrier conductiv- 
ity at temperature T = we have 



gE F r(E F ) 

2h 



(9) 



where Ep is the Fermi energy, g = 4 is the total degen- 
eracy of graphene, and r is the transport relaxation time 
at the Fermi energy obtained using the Born approxima- 
tion. The scattering time at T = due to the disorder 
potential created by charged impurities taking into ac- 
count the spatial correlations among impurities is given 



"(e P k) 



27T7J,- 



d 2 k' 

(2tt) 2 



eflk-k'l) 



x ff(#kk') [1 - cos6W] 5 (e^k' 



5(k - k') 

e P k) (10) 



where V(q) = 2ire 2 / nqe~ qd is the Fourier transformation 
of the 2D Coulomb potential created by a single charged 
impurity in an effective background dielectric constant 
K, e(q) is the static dielectric function, e s k = shvpk is 
the carrier energy for the pseudospin state "s", Vp is 
graphene Fermi velocity, k is the 2D wave vector, Ouw 
is the scattering angle between in- and out- wave vectors 
k and k', <7(#kk') = [1 + cos6*kk'] /2 is a wave function 
form-factor associated with the chiral nature of MLG 
(and is determined by its band structure). The two 
dimensional static dielectric function e(q) is calculated 
within the random phase approximation (RPA)^i, and 
given by 



Akpr s 

if q < 2kp 

— if q > 2k F 
2 



(11) 



After simplifying Eq. 1101 the relaxation time in the 
presence of correlated disorder is given by: 



n 

T 



irriihvp \ 2 f d0 (l — 



Ak f 



cos 



(sin § + 2r s ) 



fS(2fcjrsin-), (12) 



where kp is the Fermi wavevector (kp = Ep/(hvp)), 



and r s is the graphene fine structure constant (r s = 
e 2 /(hvpn) ~ 0.8 for graphene on a Si0 2 substrate). For 
uncorrclated random impurity scattering (i.e., ro = 0, 
g(r) = 1, and S(q) = 1) we recover the standard formula 
for Boltzmann conductivity by screened random charged 
impurity centers^— , where the conductivity is a linear 
function of carrier density. 

By approximating the structure factor S(2kp sin 0/2) 
that appears in (fT2"j) by a Taylor expansion around 
kp sin 9/2 = it is possible to obtain an analytical ex- 
pression for u(n) that allows us to gain some insight on 
how the spatial correlation among charged impurities af- 
fect the conductivity in MLG. Expanding the first kind 
of Bessel function J\ (x) in Eq. [5] around x ~ to the 
third order 



X 

16' 



(13) 



from Eq. (fT2j) we obtain: 



h Aimihvp 



Gi(r s )(l-Wo)+ft(r 8 )- 



H">F' 



(14) 

where the dimensionless functions G\{x) and G<z(x) are 
given by^ 3 - 



Gi{x) = - + 6x - 6ttx 2 + Ax{6x 2 - l)g(x), 



where 



4x 

— +37ra: 2 +40z 3 [l 
o 
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if X < j, 



if x > 



Using Eq. ©, 
find: 



where 



' sech 1 (2x) 
vi - 4,x 2 

sec~ 1 (2x) 
I V4x 2 - 1 

and recalling that kp 
An 

1 — a + Ba 2 n/n,i ' 



i)<?(z)L 

(15) 



(16) 



7rn, we 



(17) 
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D 



h 2n I r 2 Gi(r s ) 



H> 



(18) 



G 2 (r s ) 
2Gi(r a )' 



Note a < 1 in our model because the correlation length 
can not exceed the average impurity distance, i.e., r < 
Ti = (TOi) -1 / 2 . Eq. (|17|) indicates that at low carrier 
densities the conductivity increases linearly with n at a 
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rate that increases with ro 



whereas at large carrier densities the dependence of a on 
n becomes sublinear: 

a(n) ~ 1 - ^, (20) 
n 

where n c = (1 — a)n,/ (Ba 2 ) ~ O(l/njro). Note that 
the above equation is valid for y / 7mr <C 1, where we 
expand the structure factor as a power series of y'roro. 
The crossover density n c , where the sublinearity (n > 
n c ) manifests itself, increases strongly with decreasing 
r . This generally implies that the higher mobility an- 
nealed samples should manifest stronger nonlinearity in 
o~(n), since annealing leads to stronger impurity corre- 
lations (and hence larger ro). This behavior has been 
observed recently in experiments in which the correla- 
tion among charged impurities was controlled via ther- 
mal annealing^. Contrary to the standard-model with 
no spatial correlation among charged impurities in which 
the resistivity increases linearly in m, Eq. (|17p indicates 
that the resistivity could decrease with increasing impu- 
rity density if there are sufficient inter-impurity correla- 
tions. This is due to the fact that, for fixed ro, higher 
density of impurities are more correlated causing S(q) to 
be more strongly suppressed at low q as shown in Fig. Q] 
and [21 In the extreme case, i.e., ro = ao and n = ro, 
the charged impurity distribution would be strongly cor- 
related, indeed perfectly periodic, and the resistance, ne- 
glecting other scattering sources, would be zero. From 
Eq. (1171) we find that the resistivity reaches a maximum 
when the condition 

n/r = ^(l-TrSnr^). (21) 

is satisfied. Equation (j2"T|) can be used as a guide to im- 
prove the mobility of graphene samples in which charged 
impurities are the dominant source of disorder. 

Figs. |3l(a) and (b) present the results for <j(n) obtained 
integrating numerically the r.h.s. of Eq. (|T2"j) and keep- 
ing the full momentum dependence of the structure fac- 
tor. The solid lines show the results obtained using the 
5(q) given by the continuum model, Eq. (jSJ), the sym- 
bols show the results obtained using the 5(q) obtained 
via Monte Carlo simulations. The comparison between 
the two results shows that the analytic continuum corre- 
lation model is qualitatively and quantitatively reliable. 
It is clear that, for the same value of ro, the dirtier 
(cleaner) system shows stronger nonlinearity (linearity) 
in a fixed density range consistent with the experimental 
observations^ since the correlation effects are stronger 
for larger values of n,. 

Fig. HJa) presents that the resistivity p = 1/a in mono- 
layer graphene as a function of impurity density rii with 
correlation length ro = 5oq for different values of carrier 




FIG. 3: Calculated ff{n) in monolayer graphene with S(q) 
obtained from the Monte Carlo simulations, symbols, and 
S(q) given by Eq. solid lines for (a) m = 0.95 x 10 12 
cm -2 , and (b) m = 4.8 x 10 12 cm -2 . The different lines 
correspond to different values of ro, from top to bottom 
ro = 10ao, 8ao, 7ao, 5ao, in (a) and ro = 5ao, 4ao, 3ao, 
in (b). 




FIG. 4: (a) Calculated resistivity p in monolayer graphene as 
a function of impurity density rii for different carrier densities 
with ro = 5<io. (b) The relationship between ri/ro and y^nro 
in monolayer graphene, where the conductivity is minimum. 
The dashed line is obtained using Eq. 1211 

density. It is clear that the impurity correlations cause a 
highly nonlinear resistivity as a function of impurity den- 
sity and that this nonlinearity in p(rii) is much stronger 
for lower carrier density. In Fig. Efb) we show the value 
of the ratio r^ /ro for which p is maximum as a function of 
y/nro The analytical expression of Eq. [5l]is in very good 
agreement with the result obtained numerically using the 
full momentum dependence of S'(q). 



B. Low density: Effective medium theory 

Due to the gapless nature of the band structure, the 
presence of charged impurities induce strong carrier den- 
sity inhomogeneities in MLG and BLG. Around the Dirac 
point, the 2D graphene layer becomes a spatially inho- 
mogencous semi-metal with electron-hole puddles ran- 
domly located in the system. To characterize these in- 
homogeneities we use the Thomas-Fcrmi-Dirac (TFD) 
theory^. Ref. [H] has shown that the TFD theory cou- 
pled with the Boltzmann transport theory provides an ex- 
cellent description of the minimum conductivity around 
the Dirac point with randomly distributed Coulomb im- 
purities. We further improve this technique to calcu- 
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FIG. 5: (color online) The carrier density in monolayer 
graphene for a single disorder realization obtained from the 
TFD theory (a) for the uncorrelated case and (b) ro = 10 ao 
with m = 0.95 x 10 12 cm -2 . Carrier probability distribu- 
tion function P(n) are shown in (c), (d), (e) for (n) = 0, 
1.78, 7.7 x 10 12 cm -2 , respectively. In (f) the ratio n ruls /rii 
is shown as a function of ro/r 4 for n, = 0.95 x 10 12 cm" 2 , 
solid lines, and m = 4.8 x 10 12 cm -2 , dashed lines. We use solid lines, and m 



FIG. 6: (a) and (b) show the results for <j((n}) in monolayer 
graphene obtained from the EMT for m = 0.95 x 10 12 cm" 2 
and Hi — 4.8 x 10 12 cm -2 respectively. The different lines 
correspond to different values of ro, from top to bottom ro = 
10ao, 8ao, 7ao, 5ao, in (a) and ro = 5ao, 4ao, 3ao, in (b). 
(c) and (d) show the value of a m i„ in monolayer graphene as 
a function of ro jr{ . 



sity, (n) , and two different values of the impurity density, 
n, = 0.95 x 10 12 cm" 2 ("low impurity density") for the 



x 



10 cm ("high impurity 



(ra) = 7.7, 3.14, 0.94, x 10 cm for the solid lines (from density") for the dashed lines, 
top to bottom) and (n) = 8.34, 4.10, 1.7, x 10 12 cm" 2 for 
the dashed lines. 



late the density landscape and the minimum conductiv- 
ity of monolayer graphene in the presence of correlated 
charged impurities. To model the disorder, we have as- 
sumed that the impurities arc placed in a 2D plane at a 
distance d — 1 nm from the graphene layer. Fig. [SJa), 
(b) show the carrier density profile for a single disorder 
realization for the uncorrelated case and correlated case 
(ro = 10 ao) for = 0.95 x 10 12 cm" 2 . We can see that 
in the correlated case the amplitude of the density fluctu- 
ations is much smaller than in the uncorrelated case. The 
TFD approach is very efficient and allows the calculation 
of disorder averaged quantities such as the density root 
mean square, n rms , and the density probability distribu- 
tion P{n). Figs. EJc), (d), (e) show P(n) at the CNP, 
and away from the Dirac point (n* = 0.95 x 10 12 cm" 2 ). 
In each figure both the results for the uncorrelated case 
and the one for the correlated case are shown. P(n) for 
the correlated case is in general narrower than P(n) for 
the uncorrelated case resulting in smaller values of n rms 
as shown in Fig. [5£f) in which n lms /ni as a function of 
ro/ri is plotted for different values of the average den- 



To describe the transport properties close to the CNP 
and take into account the strong disorder-induced carrier 
density inhomogeneities we use the effective medium the- 
ory (EMT), where the conductivity is found by solving 
the following integral equatio n 2 ' 26 ' 45 " —: 



dn 



<r(n) - a emt 
V(rt) + a emt 



P{n) = 



(22) 



where u(n) is the local Boltzmann conductivity obtained 
in Section UlL&l Fig. Da) and (b) show the EMT results 
for a(n). The EMT results give similar behavior of a(n) 
at high carrier density as shown in Fig. [3j where the 
density fluctuations are strongly suppressed. However, 
close to the Dirac point, the graphene conductivity ob- 
tained using TFD-EMT approach is approximately a con- 
stant, with this constant minimum conductivity plateau 
strongly depending on the correlation length ro. Fig.[6][c) 
and (d) show the dependence of avnin on the size of the 
correlation length ro. cr m .;„ increases slowly with ro for 
r o/ r i < 0.5, but quite rapidly for ro/ri > 0.5. The re- 
sults in Fig. [HJc) and (d) are in qualitative agreement 
with the scaling of a m in with temperature, proportional 
to ro, observed in experiments^ 4 -. 
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IV. BILAYER GRAPHENE CONDUCTIVITY 

In this section we extend the theory presented in 
the previous section for monolayer graphene to bilayer 
graphenc. MLG. The most important difference between 
MLG and BLG comes from the fact that, in BLG, at low 
energies, the band dispersion is approximately parabolic 
with effective mass m ~ 0.033m e (m e being the bare 
electron mass)^ rather than linear as in MLG. As a con- 
sequence in BLG the scaling of the conductivity with 
doping, at high density, differs from the one in MLG. We 
restrict ourselves to the case in which no perpendicular 
electric field is present so that no gap is present between 
the conduction and the valence band^— . 

To characterize the spatial correlation among charged 
impurities we use the same model that we used for MLG. 



stant and is larger than 2kp for carrier density n < 
8 x 10 12 cm~ 2 . The relaxation time including correlated 
disorder is then simplified as: 



h riiTrh 2 qQ 



dx 



1 



x + qa 



2 x 2 (\-2x 2 ^ 2 



-S{2k F x) 



(27) 

where qo = qr F / '(2k F ). To incorporate analytically the 
correlation effects of charged impurities, we again expand 
S(x) around x ~ 0: 

S(2k F x) ~ 1 - a + - —a 2 x 2 - — -^a 3 x 4 (28) 
2 m 12 nf 



Combining Eqs. ([23]). ([27). and ([28]) we obtain for <r(n) 
at T = in the presence of correlated disorder 



A. High density: Boltzmann transport theory 

Within the two-band approximation, the BLG conduc- 
tivity at zero temperature T = is given by: 



(23) 



where r is the relaxation time in BLG for the case in 
which the charged impurities are spatially correlated, r is 
given by Eq. [10] with e s k = sh 2 k 2 /2m for the pseudo-spin 
state "s", e(|k— k'hthe static dielectric screening func- 
tion of BLG Ref. [56!, and ff(#kk') = [1 + cos26» kk ,] /2 
the chiral factor for states on the lowest energy bands of 
BLG. 

The full static dielectric constant of gapless BLG at 
T = is given by^ 



e(q) =[l + V(q)U(q)}- 1 

= [1 + V(q)D [g(q) - f(q)0(q - 2k F )]}~ 1 

where H(q) is the BLG static polarizability, Do ■ 
the density of states, and 



(24) 
2m 



9k 2 

m = 



2k 2 



q -^q 2 -Ak 2 F 



In 



q - vV 



1 



9{q) = W F ^ 



44 



In 



q + ^Jq 2 - 4fc 



«74 



2k p 



(25) 

To make analytical progress, we calculate the density- 
dependent conductivity using the dielectric function of 
BLG within the Thomas-Fermi approximation: 



e(q) 



1 



qrF 

q 



(26) 



e 2 2n 



1 



h n 

where 

GiM 
G 2 (qo) 
G 3 (q Q ) 



{l-a)G 1 {q Q ] + ^-a 2 G 2 [q } 



l 3 G 3 [qo] 



(29) 



9o/o 



x 2 (1 - 2x 2 



(x 


f qa 

1 


2 VT~- 

x A {\- 


- X 2 

2x 2 f 


(x 


f qo 
1 


2 VT- 

x 6 (l- 


- X 2 

2x 2 ) 2 



(X + qa ) 



dx 



dx 



dx 



(30) 



For each value of tq and carrier density n, the resistivity 
of BLG for correlated disorder is also not a linear function 
of impurity density, and its behavior is close to that in 
MLG. The maximum resistivity of BLG is found to be at 



n/r = ^2(1 - irB B 7mr 2 - C B TT 2 n 2 r^). 



(31) 



with B B G 2 [qa]/{2G 1 [q \) and C B 

— G3[(7o]/(12Gi[<7o])> which are functions weakly de- 
pending on carrier density n. 

It is straightforward to calculate the asymptotic den- 
sity dependence of BLG conductivity from the above for- 
mula and we will discuss a(n) in the strong (qo 1) and 
weak go ^ 1 screening limits separately. 

In the strong screening limit qo 3> 1, Gi[so] — tt/8, 
G 2 [q } ^ 7tt/64 and G 3 [q } ^ 13tt/128. For randomly 
distributed charged impurity, we can express the conduc- 
tivity as a linear function of carrier density u{n) ~ 
In the presence of correlated charged impurity we find: 



a(n) = 



A B n 



1-a 



,_7n_ 
16n,- 



13n 2 
192n? 



(32) 



where qrF 



4me 
ih 2 



1.0 x lO 9 !!!" 1 for bilayer graphene where a = nmr 2 , and A B 



In the strong 



on Si02 substrate, which is a density independent con- screening limit qo 1 =>■ n <C from (|32[) we ob- 



FIG. 7: Calculated o{n) in bilayer graphene with S(q) ob- 
tained from the Monte Carlo simulations (symbols) and S(q) 
given by Eq. (8| (solid lines) for two different impurity densi- 
ties (a) m = 0.95 x I0 12 cm -2 and (b) m = 4.8 x fO 12 cm" 2 . 
The different lines correspond to different values of ro. In (a) 
we use ro = 10ao, 8eto, 7eto, 5ao, (from top to bottom), 
and in (b) ro = 5ao, 4ao, 3ao, (from top to bottom). 



tain (j(n) ~ AbTi/(1 — a). With the increase of carrier 
density, the calculated conductivity in BLG also shows 
the sublinear behavior as in MLG due to the third and 
fourth terms in the denominator of Eq. 1321 

In the weak screening limit, qo <ti 1, we have Gi[go] — 
7rg§/4, G 2 [qo] and G 3 [q ] ^ 7nq$/64. The con- 

ductivity of BLG in the limit qo -C 1 is a quadratic func- 
tion of carrier density for randomly distributed Coulomb 
disorder: 



a(n) 



32n 2 



(33) 



For the correlated disorder, the calculated conductivity 
of BLG shows the sub-quadratic behavior: 



a(n) 



A b n 2 



4m 



In' 
192n 2 



(34) 



32 



with A b = — 

h niq^. F 

In Figs. [TJa) and (b), we show the a(n) within Boltz- 
mann transport theory obtained numerically taking into 
account the screening via the static dielectric function 
given by Eq. [2D We show the results for several dif- 
ferent correlation lengths ro and two different charged 
impurity densities, (a) rij = 0.95 x 10 12 cm~ 2 and (b) 
m = 4.8 x 10 12 cm~ 2 . From Figs. EJa), (b) we see that 
the conductivity increases with ro as in MLG. However 
the details of the scaling of a with doping differ between 
MLG and BLG. In BLG a(n) w n a where 1< a < 2 also 
depends on n. The effect of spatial correlations among 
impurities in BLG is to increase a at low densities and 
reduce it at high densities. 

In Fig. |U[a), we present the resistivity of BLG as a 
function of impurity density for various carrier density 
with ro = 5ao- The spatial correlation of charged im- 
purity leads to a highly non-linear function of p(rii) as 
in MLG. We also present the relation between r^/rp and 



FIG. 8: (a) The resistivity p in bilayer graphene is shown as 
a function of impurity density rii for different carrier densities 
with ro = 5ao. (b) The relationship between ri/ro and y/nro 
in bilayer graphene, where the conductivity is minimum. The 
dashed lines are obtained using Eq. 1311 



i/nro where the maximum resistivity of BLG occurs in 
Fig. EJb). The results are quite close to those of MLG 
shown in Fig. @] 



B. Low density: Effective medium theory 

As in MLG, also in BLG, because of the gapless na- 
ture of the dispersion the presence of charged impuri- 
ties induces large carrier density fluctuation s 5 ^ 57 " — that 
strongly affect the transport properties of BLG. 

Fig. [§ta) shows the calculated density landscape for 
BLG for a single disorder realization, and Fig.^a) a com- 
parison of the probability distribution function P(n) for 
BLG and MLG 57 , Within the Thomas-Fermi approxima- 
tion, approximating the low energy bands as parabolic, 
in BLG, with no spatial correlation between charged im- 
purities, P(n) is a Gaussian whose root mean square is 
independent of the doping and is given by the following 
equation 5 - 5 -: 



-l 1/2 



-f(d/r sc ) 

IT 



(35) 



where f(d/r sc ) = e 2d / r -(l + 2d/r sc )r(0, 2d/r sc ) - 1 is a 
dimensionless function, r sc = [(2e 2 m*)/(/tft 2 )] -1 rj 2 nm 
is the screening length, and T(a,x) is the incomplete 
gamma function. For small d/r sc , f = — 1 — 7 — 
log(2d/r sc ) + 0(d/r sc ) (where 7 = 0.577216 is the Eu- 
ler constant), whereas for d r sc f = l/(2d/r sc ) 2 + 
0((d/r sc )-' 3 ). As for MLG, also for BLG we find that 
the presence of spatial correlations among impurities has 
only a minor quantitative effect on P(n). For this rea- 
son, and the fact that with no correlation between the 
impurities, P(n) has a particularly simple analytical ex- 
pression, for BLG we neglect the effect of impurity spatial 
correlations on P(n). 

As in MLG the effect of the strong carrier density inho- 
mogencities on transport can be effectively taken into ac- 
count using the effective medium theory. Using Eq. (|22[) , 
cr(n) given by the Boltzmann theory, and P(n) as de- 
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FIG. 9: (Color online), (a) ra(r) of BLG at the CNP for a 
single disorder realization with n, = 10 11 cm -2 and d = 1 nm. 
(b) Disorder averaged P(n), at the CNP for BLG (MLG) 
red (blue) for m = lO^cm -2 and d = 1 nm. For MLG 
P(n = 0) ~ 0.1, out of scale. The corresponding n rms is 
5.5 x lO^cm -2 for BLG and 1.2 x lO^cm" 2 for MLG. 




FIG. 10: (a) BLG conductivity as a function of n obtained 
using the EMT for n t = 4.8xl0 12 cm" 2 for r = (4,3,2, 1,0) x 
ao from top to bottom, (b) BLG <r m m as a function of ro/n 
for m = 4.8 x 10 12 cm" 2 . 

scribed in the previous paragraph, the effective conduc- 
tivity a emt for BLG can be calculated taking into ac- 
count the presence of strong carrier density fluctuations. 
Fig. rTTJT a) shows the scaling of a with doping obtained us- 
ing the EMT for several values of tq and n, = 4.8 x 10 12 
cm . Taking account of the carrier density inhomo- 
geneities that dominates close to the charge neutrality 
point, the EMT returns a non-zero value of the conduc- 
tivity <7 m i n for zero average density, a value that depends 
on the impurity density and their spatial correlations. In 
particular, as shown in Fig- HUf b). in analogy to the MLG 
case cr m in grows with ro- 



V. DISCUSSION OF EXPERIMENTS 

Although the sublinearity of o~(n) can be explained by 
including both long- and short-range scatterers (or reso- 
nant scattcrcrs) in the Boltzmann transport theory^, it 
can not explain the observed enhancement of conductiv- 
ity with increasing annealing temperatures as observed in 
Ref. [H[ ■ Annealing leads to stronger correlations among 
the impurities since the impurities can move around to 
equilibrium sites. Our results show that by increasing ro, 
at low densities, both the conductivity and the mobility 
of MLG and BLG increase. Moreover, our results for 



MLG 32 show that as ro increases the crossover density 
at which a(n) from linear becomes sublinear decreases. 
All these features have been observed experimentally for 
MLG^. In addition, our transport theory based on the 
correlated impurity model also gives a possible expla- 
nation for the observed strong nonlinear o~(n) in sus- 
pended graphen o 21 ' 22 where the thermal/current anneal- 
ing is used routinely No experiment has so far directly 
studied the effect of increasing the spatial correlations 
among charged impurities in BLG and tested our predic- 
tions for BLG. 

Although we have used a minimal model for impurity 
correlations, using a single correlation length parameter 
ro, which captures the essential physics of correlated im- 
purity scattering, it should be straightforward to improve 
the model with more sophisticated correlation models 
if experimental information on impurity correlations be- 
comes available^. Intentional control of spatial charged 
impurity distributions or by rapid thermal annealing and 
quenching, should be a powerful tool to further increase 
mobility in monolayer and bilayer graphene devices^. 



VI. CONCLUSIONS 

In summary, we provide a novel physically motivated 
explanation for the observed sublinear scaling of the 
graphene conductivity with density at high dopings by 
showing that the inclusion of spatial correlations among 
the charged impurity locations leads to a significant sub- 
linear density dependence in the conductivity of MLG 
in contrast to the strictly lincar-in-density graphene con- 
ductivity for uncorrelated random charged impurity scat- 
tering. We also show that the spatial correlation of 
charged impurity will also enhance the mobility of BLG. 
The great merit of our theory is that it eliminates the 
need for an ad hoc zero-range defect scattering mecha- 
nism which has always been used in the standard model 
of graphene transport in order to phenomenologically ex- 
plain the high-density sublinear behavior o~(n) of MLG. 
Even though the short-range disorder is not needed to ex- 
plain the sublinear behavior of er(n) in our model we do 
not exclude the possibility of short range disorder scat- 
tering in real MLG samples, which would just add as an- 
other resistive channel with constant resistivity. Our the- 
oretical results are confirmed qualitatively by the experi- 
mental measurements presented in Ref. [44| in which the 
spatial correlations among charged impurities were modi- 
fied via thermal annealing with no change of the impurity 
density. Our results, combined with the experimental 
observation of Ref. [44[ , demonstrate that in monolayer 
and bilayer graphene samples in which charged impurities 
are the dominant source of scattering the mobility can be 
greatly enhanced by thermal/current annealing processes 
that increase the spatial correlations among the impuri- 
ties. 
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